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Abstract 



We calculate the two flavor equation of state for QCD on lattices with lattice spacing 
a = (6T)~^ and find that cutoff effects are substantially reduced compared to an earlier 
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study using a = (4T)^^. However, it is likely that significant cutoff effects remain. We 
fit the lattice data to expected forms of the free energy density for a second order phase 
transition at zero-quark- mass, which allows us to extrapolate the equation of state to 
niq = and to extract the speed of sound. We find that the equation of state depends 
weakly on the quark mass for small quark mass. 



1 Introduction 



It is generally believed that at high temperatures QCD undergoes either a phase transi- 
tion or a fairly sharp crossover into a regime where hadrons "dissolve" into a quark-gluon 
plasma. Testing this scenario is a major goal of current and planned experiments in heavy 
ion collisions. Although the fireball created in such a collision is at best in quasi-equilibrium, 
knowledge of the equilibrium equation of state for QCD is nonetheless very useful for con- 
straining the parameters of models of the quark-gluon plasma|jl|. For this reason we have 
been carrying out a program of lattice simulations to determine this equation of state, or 
energy and pressure as a function of temperature. Our calculations include two degenerate 
quark flavors. However, the strange quark is neglected, and we work at zero net baryon 
density. 

While the approximations made in a lattice simulation are controllable in principle, with 
presently available computing power these approximations are severe. In particular, effects 
of the nonzero lattice spacing are large. For example, there are big differences between the 
continuum Stefan-Boltzmann law, which is presumably the limit of the QCD energy at very 
high temperatures, and the lattice version, obtained by summing over the Fourier modes 
with the free particle action on the lattice. Also, with the Kogut-Susskind quarks that 
are usually used in high temperature QCD simulations, flavor symmetry (isospin symmetry, 
more or less) is badly broken. Instead of having N'j — 1 light pseudo-scalar particles at low 
temperature, only one pion is an exact Goldstone boson corresponding to a symmetry not 
broken by lattice artifacts. The breaking of flavor symmetry is expected to be proportional 
to (a is the lattice spacing). It can be reduced by modifying the action]^ or by simply 
decreasing the lattice spacing. 

In this work we report on an extension of our equation of state calculations to six time 
slices 0], which is a lattice spacing a = 1/(6T) instead of the a = 1/(4T) used in our earlier 
work [Q. In addition to decreasing the lattice spacing, we improve the extrapolation of our 
results to smaller quark mass by fitting our free energy to a form with either the theoretically 
predicted 0(4) critical behavior or with the mean field behavior that is expected when we 
are not very close to the critical point. 
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2 Theory 



The methods for computing the energy and pressure are standard]^, ^, 0, and we have 
discussed them in our earher paper [^]. Here we summarize the equations necessary to make 
this paper self contained. 

The energy, pressure and interaction measure are defined by 

d\ogZ 

d{l/T) 
p d log Z 

f ^ dV 

-Td\og{Z) 



E-3P 



V d\og{a) 



The temperature and volume are determined by the lattice spacing a and the space and time 
dimensions of the lattice, Ng and Nt. 

V = N^a\ 

l/T = Nta , (2) 

We use the 1x1 plaquette (Wilson) action for the gauge fields, and the conventional Kogut- 
Susskind quark action with two flavors of quarks. 

Z = J [rff/(„,^)] exp{ {6/g')Sg + {uf/A) Tr log[am, + ^] } , (3) 
Sg = \^e TTUuin,f,,u) , (4) 

n,fj,<v 

In a lattice simulation we compute derivatives of the partition function. From these we 
can either explicitly or implicitly construct the partition function and from that the energy 
and pressure. In this work we make use of these quantities: 

(□) = ^ ^^"g^^^ (5) 
^ ' 2N^Nt dQ/g^ ^ ^ 

^^^^ N^Nt darug ^ ' 

W)-C)(^^)-^g|^ (8) 
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The other second derivative of log(Z), dHog{Z) / d^aniq)'^ , involves a disconnected piece which 
is not calculated here. 

For large systems the free energy is proportional to the volume, and the pressure becomes 
just 

P=^ log(Z) (9) 

Then derivatives of the free energy are just derivatives of the pressure, and the pressure can 
be reconstructed by integrating the free energy 

^(6//,am,) = 2N^N, f''\lu{Q/g'\am,)) - (0(6/9", am,)) ]d{6/g") , (10) 

1 J cold ^ ' ^ ' sym 

or 

^{Q/g^ am,) = N^Nt / [(#(6//, m'^a)) - U^{Q/g\ m'^a)) ]d{m'^a) , (11) 

1 v/cold ^ ' ^ ' sym 

where the "symmetric" quantities subtract the divergent zero-temperature pressure. The 
interaction measure can be found from simulations at a single value of (6/(7^, amg) and the 
beta function, which tells how these lattice couplings must be changed to change the lattice 
spacing. 

'-^ = - 2iV|A',«|C> - - N':nS^[{m,) - {M,) ] (12) 

1 a log a ^ a log a ^ ' ^ ' sym 

In reference!^ we obtained the nonperturbative beta function for couplings associated with 
the Nt = 4 and 6 crossovers. To briefly summarize, zero temperature spectrum data from the 
literature are combined in a fit which gives am^ and amp as smooth functions of Q/g"^ and 
am,. The p mass is used to define the lattice spacing. In particular, we somewhat arbitrarily 
set the p mass to 770 MeV at all light quark masses. Lines of constant physics {i.e. , m-^/mp = 
constant), along which the lattice spacing varies, are determined in the bare coupling space. 
The two components of the beta function, d{6/ g"^)/ (loga) and d{amq)/d (log a), tell how the 
input parameters change along these lines of constant physics. 



3 Simulations 



For the equation of state, we have carried out simulations using the hybrid molecular dy- 
namics R algorithmp]]. The calculation requires both asymmetric (12'^ x 6) and symmetric 
(12^) lattices as mentioned above. These are commonly referred to as hot and cold lattices, 
respectively, and we will use this terminology. This is somewhat misleading since the system 
may be in the cold phase, i.e. below the crossover, even on the hot lattices if the coupling 
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is small enough. Each hot (cold) simulation is at least 1800 (800) time units long after 
equilibration. As a rule, fewer trajectories were required for the "cold" lattices to achieve 
the same level of statistical accuracy because the lattice volume was twice as large. On the 
hot lattices near the crossover region the simulations were extended to more than 3000 units 
to overcome the larger fluctuations associated with the transition. 

The bare coupling and quark mass phase diagram corresponding to our simulations is 
shown in Fig. 0. The vertical bars denote the approximate location of the finite temperature 
crossover for airtq = 0.0125 and 0.025. Note, while increasing the coupling Q/g"^ is analo- 
gous to increasing the temperature, lowering the bare quark mass aruq also increases the 
temperature. (This statement depends on the physical quantity used to define the lattice 
spacing. Again, we use mp.) Of course there are lines in this bare coupling phase diagram 
on which m.,^/mp is fixed and only the physical temperature of the system varies. In par- 
ticular, the bottom of the graph, where m^/mp = 0, is the line arug = and corresponds 
approximately to the real world. In our simulations the gauge coupling takes the values 
5.35 < Q/g"^ < 5.6 for quark mass arUq = 0.0125 and 5.37 < Q/g"^ < 5.53 for aniq = 0.025. 
Then along the line Q/g'^ = 5.45, the quark mass varies in the range 0.01 < arrig < 0.1 
and at 6/g'^ = 5.53, 0.0125 < arriq < 0.2. This range of couplings and masses corresponds 
roughly to physical temperatures 125 < T < 250 MeV (based on the p mass[0) and mass 
ratios 0.3 < m^/mp < 0.7. In units of the pseudo-critical temperature Tc, m^/Tc = 1.94 
and 2.69 for arUq = 0.0125 and 0.025, respectively. Past lattice simulations indicate that 
Tc/rrip ^ 0.20, or itLt^/Tc ~ 0.9 in the continuum limit. Thus our simulations correspond to 
rather heavy pions. 

The R algorithm introduces lowest order errors in observables that are proportional to 
At^, where At is the step size used to numerically integrate the gauge field equations of 
motion through simulation time^j. The errors are in general different on hot and cold 
lattices, thus multiple simulations at each value of Q/ g'^ and airiq are required to extrapolate 
observables to At = 00]. Of course, this greatly increases the computational cost of the 
calculation since at each time step the force term due to the quarks in the gauge field 
equations of motion requires the inverse of the quark matrix, and as it happens. At must be 
taken relatively small (~ auiq) to be in the regime where the lowest order error dominates. 
Thus many inversions of the quark matrix are required for each simulation. In Fig. ^ we 
show example results for the plaquette ((□)) at airiq = 0.0125. Evidently, on the cold lattices 
and for lower values of Q/g"^, the effects are worse. The step size errors are particularly 
troublesome for the plaquette since the difference of the plaquette on hot and cold lattices is 
quite small (even in the hot phase), as we will see in the next section. At low temperature 
both (ipip) and (□) approach their symmetric values, making extraction of thermodynamic 
quantities in the hadronic phase very difficult. 
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4 Results and Analysis 



The expectation values of (□) and {^ip) are shown in Figs. ^ and ^ These values reflect the 
step size extrapolations. In each figure, a small but noticeable inflection point is observable in 
the expectation values on the hot lattices at 6/(?^ ~ 5.415(amg = 0.0125) and 5.445(0771^ = 
0.025). These couplings correspond to the pseudo-critical temperatures of the crossover. 
The expectation values are smooth on the cold lattices. Note that the values on the hot 
lattices approach their cold values at small Q/g^ and begin to separate as Q/g^ (and thus 
the temperature) is increased. The area between the curves in Fig. |^ yields p/T^- At large 
6/(7^, the plaquette expectation values again approach each other since p/T'^ —>■ constant as 
T —>■ oo, as expected from asymptotic freedom of QCD. The qualitative behavior of (ipijj) 
is consistent with its expected role as an order parameter for a spontaneously broken chiral 
symmetry. A quantitative analysis and extrapolation to zero-quark-mass is given later. 

Next we turn to a discussion of the pressure. The integration of (V^V') with respect to arug 
at Q/g'^ = 5.45 and 5.53 yields the pressure as a function of aruq at fixed 6/g'^, shown in Fig. |[ 
For large arrig the system is in the cold phase, and as arUg decreases, the pressure smoothly 
increases. To extrapolate the pressure from our smallest quark mass to zero-quark-mass 
requires a corresponding extrapolation of the cold and hot contributions to the integrand 
in ([TT|). For the hot lattices at both values of Q/g^ the extrapolation takes place in the hot 
phase, so for this contribution we assume {i'i') = at amg = 0. For the cold lattices we 
extrapolate using the fit summarized in Table |l|. A simple linear extrapolation of the hot 
lattice data for the smallest two quark masses gives a result within 1.5a of zero at Q/g^ = 5.45 
and 3.0(7 of zero at Q/g"^ = 5.53, perhaps due to curvature effects or underestimated errors. 
In any event forcing a zero intercept has a negligible effect on the pressure extrapolation. 

As mentioned above, the pressure as a function of Q/g"^ at fixed quark mass is obtained 
by integrating (□) with respect to Q/g'^. In Fig. ^ we show the results for arUg = 0.025 
and 0.0125. Again, the pressure rises smoothly through the crossover. The curves are 
similar except for an overall shift in G/g"^, the arrig = 0.0125 curve beginning its rise sooner 
since the smaller quark mass corresponds to higher temperature. Shown where they can be 
compared are the points obtained from the the quark mass integration. The values from the 
two different approaches agree, indicating that the integration method works well for the 
volumes studied. 

The interaction measure at each point is just the sum of the T = subtracted values of 
the (□) and (ipip) weighted by the coupling and quark mass components of the (3 function, 
respectively (Eq. |12|). The results for arrig = 0.0125 and 0.025 are shown in Fig. |^. At zero 
temperature the interaction measure is zero since we have normalized e and p to be zero at 
T = 0. Through the transition, we expect I to increase rapidly if the energy density increases 
rapidly, e.g. if the quarks and gluons deconfine, since the pressure must rise smoothly (even 
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for a discontinuous transition) to maintain mechanical equilibrium of the system. This 
behavior is seen in Fig. |^. Again, at higher temperature due to asymptotic freedom, we 
expect / to decrease to zero as both the energy density and pressure asymptotically approach 
their Stefan-Boltzman values, and e = 3p, the equation of state for a relativistic free gas. 

The energy density constructed from I and p is shown in Fig. |[ Here we include 3p as 
well. We also plot e and 3p, using the values of the (□) and {ipip) at the smallest step size 
available at each point, without step size extrapolations. The difference between the two 
gives an estimate of the step size systematic error in our final results. Recall that step size 
effects for larger masses were essentially eliminated by taking At -C aniq. 

In Fig. ^ we show the equation of state as a function of the physical temperature. The 
quark mass dependence is largely removed by this partial rescaling, {6/g'^{a),amq{a)} 
{T{6/g'^,amg),amg}. Note that the physical quark mass, or more precisely mT^/mp, varies 
along each line of constant aruq while m^/T = arUgNt is held fixed. From the figure, we see 
a large increase in the energy density as T increases through the pseudo-critical temperature 
Tc ~ 140 MeV {Q/g"^ was defined above). Just below T^, e/T^ ^ 5 which is already substan- 
tial. Physically this means the hadronic phase has non-negligible energy density, except that 
it is not clear what degrees of freedom are being excited since we have already mentioned 
that m^/Tc is roughly 2 to 3. A similar situation exists for the SU(3) pure gauge case where 
the ratio of the lightest glueball state to Tc is greater than 5[0. We also note that the 
energy density of a relativistic gas of three light pions is insufficient to explain the observed 
energy density in the hadronic phase {e/T^ = vr^/SO x 3 ~ 1). Similarly, in SU(2) pure 
gauge theory it has been noted that a gas of non-interacting glueballs cannot account for 
the observed energy density below Te [p]] . 



In Fig. 1^ we also compare the Nf = 6 equation of state with an earlier result on A^^^ = 4 
lattices and with the continuum and lattice Stefan-Boltzmann laws. There is an apparent 
large finite size effect which is expected from the free lattice theory. At high temperature the 
energy density has leveled off dramatically while the pressure is still increasing at the largest 
value of T that we simulated. Because of asymptotic freedom, e and 3p should approach the 
Stefan-Boltzmann result for the corresponding value of Nt. But, from Fig. ^ the approach 
to the free result is evidently quite slow. 



4.1 Extrapolation to ruq = 

The above results pertain to unphysical values of the quark mass. Indeed, we would like to 
obtain the equation of state along a line of constant physics corresponding to the real (two 
flavor!) world. This can be done by extrapolating to the chiral limit, arrig 0. To this end 
we fit the derivatives of the free energy density to an appropriate function of Q/g'^ and arriq 



7 



and then set aruq = 0. The fit serves the dual purpose of smoothing the data and allowing 
a parameterization of the equation of state in terms of the bare quantities, from which we 
can extract, e. g., the speed of sound. 



If the QCD high temperature phase transition with two flavors is a second order transition 



at zero- quark- mass driven by the restoration of chiral symmetry! 12], then the critical part 



of the free energy should have a universal form, up to the scale of the gauge coupling and 



quark mass|13]. The free energy is the sum of an analytic piece and a scaling piece: 

f = fa{Q/g^am,) + f,{t,h) (13) 

where t = {T — Tc)/Tq and h = H / Hq. Tq and Hq are conventionally determined by requiring 
{'4)il)){t = 0,h) = h^^^ and {ijjip){t < 0,h = 0) = (—t)^. (In the language of spin models, 
(ipip) is the magnetization and (□) is the energy.) From invariance under a length rescaling 
by a factor b, the critical part of the free energy should have the property: 

f,{t,h) = h-''f,{hy%hy>^h) . (14) 

This implies that the magnetization near the critical point is determined by a universal 
scaling function, conventionally written as: 

^ fit/h'/f'') = fix). (15) 



The normalization conditions on t and h then require that /(O) = 1 and f{x) {—x)^ as 
X ^ —oo. This condition, along with the known values of i/h and ?/t[ll4|, has been used to 
compare the behavior in (ijjijj) with that expected from 0(4) symmetry [p!5|]. Here we wish 
to use this theoretical input to guide extrapolation of the free energy to physical light quark 
mass, which is essentially the same as arUq = 0. We fit to both a scaling form for 0(4) in 
three dimensions and to the form for mean field theory. We expect the mean field form to be 
a good approximation when the system is not very near the critical point, and the difference 
between these two forms gives an idea of the systematic errors in this approach. 

Since we calculate both {ipip) and (□) (and their derivatives), we would like to treat 
them equally in fitting the free energy. Therefore we use a formulation of the scaling free 
energy which handles the energy and magnetization symmetrically. This has been discussed 



in Ref. |[16||, so we just summarize it here. 



The scaling ansatz, Eq. |T4|, tells us that if we specify the singular free energy on the unit 
circle in the t, h plane, we have specified it for all t, h. Thus the scaling part of the free 
energy density can be written 

fsit,h) = bit,h)-'g{eit,h)) (16) 

where b is the solution to 

(bynf + {W'^hf = 1 (17) 
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fit Q/gl To H, Us xVdof 

0(4) 5.353 0.522 1.51 3.339A/5 + 2.30A/32 + 0.98A/?3 

+am2(0.41 + 4.3A/3-8.3A/?2) 990/91 
MF 5.381 0.546 1.22 3.221A/5 + 0.505A/52 - 0.221A/53 

+am2(0.63 + 0.7A/?-0.6A/32) 933/91 
cold 0.132A/3 + 2.02A/52 + 3.195A/33 + -1.23A/?4 130/52 

+amg{lMAP + 0.723A/?2 + 5.5A/33 - 45A/?4) 
+am^(-1.34A/? - 12A/3^ + 215A/3^ + 2.2A/?^) 

Table 1: Fit summary table, fns is the non-singular free energy up to an overall constant. 
Note A/3 = Q/g^ - 5.4. 



and 



e{t,h) =tan-^(6^^/i,6^*t) 



(18) 



Here g{9) is a universal function. In Ref. |[T6[| g{6) for 0(4) is determined approximately by 
Monte Carlo simulation of the 0(4) spin model. For this formulation it is convenient to 
modify the conventional normalizations of t and h, and use t = {T — Tc)/Tg and h = H/Hg, 
where Tg and Hg are chosen to fix 5'(vr/2) = yh/d and g'{7i) = —1. These are related to the 
conventional normalizations by Hq = H^'^^ and Tq 



From simulation of the 0(4) spin model an approximate g{6) for 0(4)[|TB] has been 
obtained. For the mean field case, the scaling function can be obtained from a numerical 
re-parameterization of the mean field magnetic equation of state [TB[ 



(19) 



Fig. ^ shows (ip'ip) calculated from fits to the mean field and 0(4) scaling functions, plus 
polynomials in arriq and 6/g^. We also include the pure polynomial fit to the cold data, and 
both hot and cold data are shown for comparison. The fits are summarized in Table P^r 
degree of freedom is poor for all of the fits. For the hot data, the mean field and 0(4) cases 
each have x^/dof ^ 1000/91 while the cold data has x^/dof ^ 130/52. If we fit to the data 
without step size extrapolations, the results improve somewhat, x^/dof ~ 650/91 for the hot 
data and 117/56 for the cold. A fit to only (iptp) over the same range gives x^/dof ~ 49/18. 
Despite the high x^, examination of Fig. ^ shows that the data are actually reproduced quite 
well by the fits. This is remarkable given the large range of arriq and Q/g"^ spanned by the 
fits. Moreover, the fits to mean field and 0(4) scaling functions work equally well; our data 
cannot distinguish between the two scaling behaviors. However, it is interesting to note that 
the respective extrapolations to rriq = are quite different, and give critical temperatures 
Tc ~ 140 and 150 MeV for 0(4) and mean field respectively. The above indicates that the 
present lattice simulations may still be too far from the scaling region and smaller quark 
masses are required to see the true scaling behavior. The correlation length in lattice units 
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as given by the inverse pion mass (the pion and sigma are degenerate at the critical point) 
of the present simulations is only 2-3 which is less than the temporal extent of the lattice, 
so true dimensional reduction, which must occur for universality arguments to hold, has not 
been achieved. 

At arriq = 0, d{amq)/d{loga) = so the interaction measure is determined solely from 
the plaquette. However, most of our information about the critical behavior comes from 
(ipip), since the contribution of the scaling part of the free energy to the plaquette is small 
compared with the analytic part. Therefore in a crude approximation our procedure is using 
information about (ipip) to help determine the free energy, which in turn yields the plaquette 
as arriq 0. 



An extrapolation of the equation of state to niq = is shown in Fig. [Ty. It is compared 
to the arUq = 0.0125 result which reproduces the data reasonably well. The appearance 
of the bump in the energy density just after the transition is probably an artifact of the 
extrapolation (at rriq = 0, the corresponding region of G/g"^ lies below the values of the 
coupling where we have done simulations). From Fig. [1^ we again see a weak dependence on 
the quark mass, which gives us some reassurance in the extrapolation. The plus and minus 



one standard deviations shown in Fig. |TO]are calculated in the following way. First, we do a 
covariant fit to (□), (ipip), and their derivatives with respect to Q/g"^ (Eqs. §-§). From the 
fit we obtain a set of parameters and the covariances of these parameters which map onto a 
multi-dimensional Gaussian probability distribution for the parameters. We then generate 
many parameter sets with this distribution and calculate the equation of state for each one. 
The standard deviation of the mean of this set is shown in the figure. 



In Fig. |Tl| we show the speed of sound squared calculated from the 0(4) fit. It rises 



rapidly through the transition region and then levels off near the free value of 1/3. This 
indicates that the system is weakly interacting in this region. This should be contrasted 
with e — 3p just after the crossover which indicates significant interaction effects. Indeed, the 
couplings in the region are of order one, and we have already seen that neither the energy 
density nor pressure is approaching its perturbative value. The low temperature part of the 
curve is probably not accurate. The derivatives of the energy density and pressure are poorly 
known in this region since the difference of (□) and (ijjip) from their cold lattice values is 
nearly zero. We have already mentioned that m.,^/Tc was rather high in our simulations, so 
it is not surprising that the expected dynamics of a dilute gas of relativistic pions is not 
observable. In that case, we expect the hadron gas below the transition to have a nonzero 
speed of sound, which then dips down at the transition. The statistical errors for the speed 
of sound were calculated in the same manner described above for the equation of state. 
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5 Conclusions 



In this work we have pushed our calculation of the equation of state for QCD including 
dynamical quarks to smaller lattice spacing. We have also developed techniques for using 
theoretical expectations for the scaling behavior to extrapolate numerical results to the 
physical quark mass. The points where we ran were chosen to explore the equation of state 
over a range of temperature rather than to work very close to the crossover. Therefore, we can 
say little about the nature of the transition or crossover from this work. Still, our results are 
consistent with the standard picture of a second order phase transition at zero-quark-mass 
and a sharp crossover for small but nonzero masses. 



We remark, however, that recent simulations with two flavors of Kogut-Susskind quarks 
on lattices with A^^ = 4 and at smaller quark masses than used in this study have revealed 
significant finite size effects that have cast doubt on earlier promising demonstrations of 



critical scaling [|T^. Other recent simulations with two flavors of Wilson quarks and an 
improved gauge action to reduce cutoff effects also show promising agreement with 0(4) 



scaling, but a thorough investigation of finite-size effects remains to be done||18||. So for 
the moment the question of the order and universality class of the transition remains open, 
and the validity of our extrapolation to zero quark mass remains to be established. At the 
very least, one may expect that the recently reported finite volume corrections lead to a 
greater sharpening of the crossover in energy density and speed of sound as the quark mass 
is decreased. We have seen that small changes in the extrapolation have a large effect on 
some but not all extrapolated values: e.g. 0(4) and mean-field extrapolations give zero- 
quark-mass critical temperatures that differ by 10 MeV. On the other hand, expressed as 
a function of temperature in physical units, the energy density and pressure away from the 
crossover show little dependence on quark mass, even in the zero-quark- mass limit. 



The methods used here require a subtraction of the zero temperature plaquette and 
{ipip)- As Nf is increased, the plaquette subtraction rapidly becomes more difficult, since 
the fractional difference in plaquette between the hot and cold lattices decreases as A^^^"^. 
This suggests that, as for many other quantities, an improved action which allows use of 
a larger lattice spacing, or smaller A^^, will be important for further progress. Results for 
pure gauge theory and for four-flavor QCD have been reported by the Bielefeld group [p^. 
Another important problem for future studies is remedying the breaking of flavor symmetry, 
so that the low temperature phase that is simulated really has three light pions. For a start 
in this direction see Ref. p|. 
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Figure 1: Phase diagram for our simulations. The vertical solid lines indicate approximate 
locations of the crossover. Crosses indicate cold lattices, octagons hot lattices. 
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Figure 2: The plaquette expectation value as a function of the step size squared, Ai^, and 
the gauge coupUng. Results are shown for arriq — 0.0125. Octagons denote hot lattices 
is given for each) and squares denote cold lattices. 
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Figure 3: The plaquette expectation value as a function of the gauge couphng at arriq = 
0.0125. Octagons denote hot lattices, squares cold lattices. The area between the curves 
gives the pressure while the difference at each point is related to the interaction measure. 
The lines are from a fit including the 0(4) singular free energy described in the text (dashed 
lines are an extrapolation to arUg = 0). 
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Figure 4: {ipip) as a function of Q/ and am^. The difference in the hot (octagons) and cold 
(squares) values is related to the interaction measure, {ipip) on the hot lattices also serves 
as an order parameter for the system. Lines depict fits to 0(4) and mean field singular 
forms of the free energy plus analytic terms to (□), (ipip), and their derivatives (Eqs. |^ 
Only points with arriq < 0.025 were included in the fits. The difference between 0(4) 
(dotted line) and mean field forms (dashed line) is not discernible at the quark masses where 
simulations were run. Extrapolations to rriq = 0, however, give different critical couplings. 
Solid lines correspond to polynomial fits to the cold data and the corresponding extrapolation 
to arriq = 0. 
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Figure 5: The pressure obtained from integrating {ipil') with respect to aniq and constant 
G/g"^. The values at zero- quark-mass (bursts) are obtained by setting {'ipip) = on the hot 
lattices and extrapolating {ipip) to auig = on the cold lattices. 
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Figure 6: The pressure obtained from integrating (□) with respect to 6/(?^ at constant aniq. 
The values from the quark mass integrations (fancy squares) are also shown for comparison. 
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Figure 7: The interaction measure, s — 3p. 
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Figure 8: The energy density constructed from I + 3p (upper two curves). Also shown is 3p 
(lower two curves) and results with no step size extrapolations (squares). The difference is 
the step size systematic error in the equation of state. The data are for aniq = 0.0125 only. 
The fancy squares denote 3 times the pressure as calculated from the integration of {ipip) 
with respect to arUg. 
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Figure 9: The equation of state along lines of constant rUq/T. The octagons denote 
aruq — 0.0125 results, squares aruq — 0.025. The diamonds indicate an earlier result on 
Nt — 4: lattices. Horizontal lines correspond to Stefan- Boltzmann laws for Nt = 4, 6, and 
the continuum. The energy density increases rapidly near the crossover while the pressure 
(lower curve for each symbol) rises smoothly. 
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Figure 10: The arriq = equation of state from a fit to the data that includes the 0(4) 
universal scaling function. Also shown is the data and the fit at aniq — 0.0125 for comparison. 
Again, there is only a weak mass dependence. The "bump" just after the transition is likely 
an artifact of the extrapolation. The solid lines correspond to the central value and a one 
standard deviation above and below this result (statistical errors only). 
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Figure 11: The speed of sound squared from a fit to the data that includes the 0(4) 
universal scaling function, rises rapidly in the crossover region and then approaches the 
free gas result, 1/3. The low temperature result is probably not accurate (see text). The 
solid lines correspond to the central value and a one standard deviation above and below 
this result (statistical errors only). 
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